Universal dynamics in the onset of a Hagen— Poiseuille flow 
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The dynamics in the onset of a Hagen-Poiseuille flow of an incompressible liquid in a channel of 
circular cross section is well-studied theoretically. We use an eigenfunction expansion in a Hilbert 
space formalism to generalize the results to channels of arbitrary cross section. We And that the 
steady state is reached after a characteristic time scale r = {A/'P)'^{l/u) where A and P are the 
cross-sectional area and perimeter, respectively, and v is the kinematic viscosity of the liquid. For 
the initial dynamics of the flow rate Q for t ^ r we flnd a universal linear dependence, Q{t) = 
Qoo {ct/C) (t/r), where Qao is the asymptotic steady-state flow rate, a is the geometrical correction 
factor, and C = /A is the compactness parameter. For the long-time dynamics Q{t) approaches 
Qao exponentially on the timescale r, but with a weakly geometry-dependent prefactor of order 
unity, determined by the lowest eigenvalue of the Helmholz equation. 

PACS numbers: 47.10.A-, 47.15.Rq, 47.27.nd, 47.27.nf, 47.61.-k, 47.85.L- 



I. INTRODUCTION 

Hagen-Poiseuille flow (or simply Poiseuille flow) is im- 
portant to a variety of applications ranging from macro- 
scopic pipes in chemical plants to the flow of blood in 
veins. However, the rapid development in the field of 
lab-on-a-chip systems during the past decade has put 
even more emphasis on pressure driven laminar flow. 
Traditionally, capillary tubes would have circular cross- 
sections, but today microfabricated channels come with a 
variety of shapes depending on the fabrication technique 
in use. The list of examples includes rectangular chan- 
nels obtained by hot embossing in polymer wafers, semi- 
circular channels in isotropically etched surfaces, trian- 
gular channels in KOH-etched silicon crystals, Gaussian- 
shaped channels in laser-ablated polymer films, and ellip- 
tic channels in stretched PDMS devices, see e.g. Ref. 01 

This development has naturally led to more empha- 
sis on theoretical studies of shape-dependence in mi- 
crofluidic channels. Recently, we therefore revisited the 
problem of Poiseuille flow and its shape dependence and 
we have also addressed mass diffusion in microchan- 
nels 1^ 1^. In the present work we combine the two 
former studies and address the dynamics caused by the 
abrupt onset of a pressure gradient at time i = in an 
incompressible liquid of viscosity ij and density p situ- 
ated in a long, straight, and rigid channel of length L 
and some constant cross-sectional shape Q. The solution 
is well-known for the case of a cylindrical channel, see e.g. 
Ref. 0, but in this paper we generalize the results to a 
cross-section of arbitrary shape. The similarity between 
mass and momentum diffusion, and the existence of a 
characteristic diffusion time-scale r^^g = {tt / 'i){A/V)'^ / D 
for mass diffusion [3] have led us to introduce the momen- 
tum diffusion time-scale r defined by 



(1) 



sions of a diffusion constant), while A and V is the area 
and perimeter of the cross section fi, respectively. In this 
paper we show that the dynamics of the flow rate Q{t) is 
universal with t as the characteristic time scale. 

As shown in Ref. the shape parameters A and V also 
play an important role in the steady-state Poiseuille flow. 
The hydraulic resistance -Rhyd can be expressed as 



U '^■^ - 73* 



(2) 



where a is a dimensionless geometrical correction factor 
and -RJ^y^ = rjL/A^ is a characteristic resistance. Re- 
markably, a is simply (linearly) related to the dimen- 
sionless compactness parameter C = V"^ / A. 

Above we have emphasized microfluidic flows because 
of the variety of shapes frequently encountered in lab-on- 
a-chip systems. However, our results are generally valid 
for laminar flows at any length scale. 



II. DIFFUSION OF MOMENTUM 

We consider a long, straight channel of length L, 
aligned with the z-axis, having a constant cross section 

with the boundary dVL in the xy plane. The fluid flow 
is driven by a pressure gradient of Vp = —{Ap/L)ej; 
which is turned on abruptly at time t — 0. We note 
that strictly speaking the pressure gradient is not es- 
tablished instantaneously, but rather on a time-scale 
set by L/c where c is the speed of sound. For typi- 
cal liquids c ^ lO'^m/s which for micro-fluidic systems 
and practical purposes makes the pressure gradient ap- 
pear almost instantaneously. From the symmetry of the 
problem it follows that the velocity held is of the form 



v(r,t) = v{rj_,t)ex where r_L 



yey 



ze,. From the 



Navier-Stokes equation it then follows that v(r±,t) is 
governed by, see e.g. Refs.|3or|^ 



where v — rj/ p is the kinematic viscosity (having dimen- 



■dtv{v^,t)-V^v{v^,t) 



Ap 



(3) 
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which is a diffusion equation for the momentum with the 
pressure drop acting as a source term on the right-hand 
side. The velocity v is subject to a no-sHp boundary 
condition on 951 and obviously v is initially zero, while it 
asymptotically approaches the steady-state velocity field 
i'oo(r±) for t — !■ oo. 

In the analysis it is natural to write the velocity as a 
difference 



v{rj_,t) = Voo{r±) - v,^{rj_,t) 



(4) 



of the asymptotic, static field Wqo, solving the Poiseuille 
problem 



-V\oo(r±) = 



Ap 



(5) 



and a time-dependent field v^^ (rj_ , t) satisfying the homo- 
geneous diffusion equation, 



-atfJri,t)-VV(r-L,i)-0. 



(6) 



From Ref. 3 it is known that rescaling the Helmholz equa- 
tion by {A/V)"^ leads to a lowest eigenvalue ai that is of 
order unity and only weakly geometry dependent. We 
therefore perform this rescaling, which naturally implies 
the time-scale r of Eq. JP) and the following form of the 
diffusion equation, 

Tdtv^{r^,t)~^v^{r^,t) = 0. (7) 

where we have introduced the rescaled Laplacian Jif, 

(8) 



We note that by the rescaling the Navier-Stokes equa- 
tion (O becomes 



r dtv — — ^ 



(9) 



where we have introduced the steady-state flow rate 
Qoo = Ap/_Rhyd and used Eq. 0. 



III. HILBERT SPACE FORMULATION 

In order to solve Eq. ^ we will take advantage of the 
Hilbert space formulation ^ , often employed in quantum 
mechanics 0- The Hilbert space of real functions f{r±) 
is defined by the inner product 



{f\g)^ [ drx /(rx)g(r^) 



(10) 



and a complete set {| of orthonormal basis func- 
tions, 



(11) 



Above, we have used the Dirac bra-ket notation and Snm 
is the Kronecker delta. We choose the eigenfunctions 
{|0n)} of the rescaled Helmholz equation (with a zero 
Dirichlct boundary condition on 9f2) as our basis func- 
tions. 



-if l^n) = a„|0„). 



(12) 



With this complete basis any function in the Hilbert 
space can be written as a linear combination of basis 
fimctions. Using the bra-ket notation Eq. Q becomes 



■dt\v)-^\v)^^\l). 



(13) 



The full solution Eq. is written as 

\v) = \v^)-\vh), (14) 
where |foo) satisfies the Poiseuille problem Eq. (O, 



-=^K) = ^|i), 



(15) 



and the homogeneous solution solves the diffusion 
problem Eq. 



{rdt~^)\v,,)^0. 
In the complete basis {|(/>n)} we have 

oo 

n=l 
oc 

\voo) = E C„|0„), 



(16) 

(17) 
(18) 



and since limj_+o = ["^cx;) we have 6„ = c„. Multi- 
plying Eq. (UHl) by (^,„| yields 



Cm = (0m |woo) = (0m l-Sf ^^\voc) 
= ^^-flm (0m 1 1}- 



(19) 



In the second-last equality we have introduced the unit 
operator 1 — and in the last equality we used 

the Hermitian property of the inverse Laplacian opera- 



tor to let ^ ^ act to the left. 



from Eq. l(T^ . while ^ acts to the right, see Eq. l(T^ . 
Substituting Eqs. (^71) and I^Hl into Eq. we finally 
obtain 



-1 -a„t/'- 



(20) 



IV. FLOW RATE 



Using the bra-ket notation, the flow rate Q(t) at any 
time is conveniently written as Q = (ij^), and thus in 
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steady state Qoo = (l|woo)- Multiplying Eq. from 
the left by (l| yields 

Q{t) = {l\v) = Qoo-^ E (lkn)('^«|l)<'e--*/^ 

n— 1 

(21) 

The factor (l|0„)((/)„|l) is recognized as the effective 
area An covered by the nth eigenfunction |0„), 

An ^ = Kl^'Or = (l|</'n)(^„|l). (22) 

The effective areas fulfil the sum-rule — A, 

seen by completeness of the basis {|0n)}- We can find 
the geometrical correction factor a from Eq. (|21|l by using 
that Q(0) = 0, 

""iS-)' '''' 

and substituting into Eq. 121|) we finally get 

Qoc \^^^ an J an 

V. SHORT-TIME DYNAMICS 



shape 


ax 


Ai/A 


a/C 


circle 


7g,i/4 ~ 1.45'' 


4/7^,1 ~ 0.69" 


2" 


quarter-circle 


1.27^= 


0.65" 


1.85" 


half-circle 


l.SS'^ 


0.64" 


1.97" 


ellipse(l:2) 


1.50" 


0.67" 


2.10" 


ellipse (1:3) 


1.54^= 


0.62" 


2.21" 


ellipse(l:4) 


1.57^= 


0.58" 


2.28" 


triangle(l:l:l) 


7rV9 ~ 1.10'* 


~ 0.61'* 


5/3* 


triangle(l:l:\/2) 


4(2+72)2 - -^-^0 


512/97r'' ~ 0.58" 


1.64" 


square 


ttVS ~ 1.23" 


64/7r* ~ 0.66" 


1.78" 


rectangle(l:2) 


57rV36 ~ 1.37" 


64/7r* ~ 0.66" 


1.94" 


rectangle(l:3) 


57rV32 ~ 1.54" 


64/7r* ~ 0.66" 


2.14" 


rectangle(l:4) 


177rVl00 ~ 1.68" 


64/7r* ~ 0.66" 


2.28" 


rect angle (licxa) 


~ 7rV4 ~ 2.47" 


64/7r* ~ 0.66" 


~3" 


pentagon 


1.30" 


0.67" 


1.84" 


hexagon 


1.34" 


0.68" 


1.88" 



TABLE I: Central parameters for the lowest eigenfunction 
for different cross sectional shapes. Note how the different 
numbers converge when going through the regular polygons 
starting from the triangle(l:l;l) through the square, the pen- 
tagon, and the hexagon to the circle. "See e.g. Ref. |g for 
the eigenmodes and eigenspectrum. Here, 70,1 denotes the 
first root of the zeroth Bessel function of the first kind. ''See 
Ref. "Data obtained by finite-element simulations. '*See 
e.g. Ref . 1^ for the eigenmodes and eigenspectrum. "See e.g. 
Ref. U for a solution of the Poisson equation. 



The short-time dynamics is found by Taylor-expanding 
Eq. |(23|) to first order, 



Q(t) 



An 



. n=l 



t aA t at 



t < r. 



(25) 

where we have used the sum-rule for An as well as 
Eq. (|23|l . The short time dynamics can also be inferred 
directly by integration of the Navier-Stokes equation 
Eq. H13|l . since at time t = we have |?;) = and con- 
sequently the vanishing of velocity gradients and viscous 
friction, .^\v) = 0. Thus we arrive at 



T dt\v) = 



-p2 



(26) 



corresponding to a constant initial acceleration through- 
out the fluid. Integration with respect to t is straight- 
forward and multiplying the resulting |u) by (l| yields 

0(t), 



Q{t) 

Qoo 



a / I N t aA t at 



t < T. 



(27) 



Thus initially, the fluid responds to the pressure gradient 
in the same way as a rigid body responds to a constant 
force. 



VI. LONG-TIME DYNAMICS 

As the flow rate increases, friction sets in, and in the 
long-time limit t » t the flow-rate saturates at the value 
Qoo where there is a balance between the pressure gra- 
dient and frictional forces. For the long-time saturation 
dynamics the lowest eigenstate plays the dominating role 
and taking only the n = 1 term in Eq. H24|l we obtain 



Q{t) 



l-e-'^l*/^ t>T/a2, 



(28) 



where we have used that the lowest eigenvalue ai is non- 
degenerate to truncate the summation. 

The time it takes to reach steady-state is denoted Too- 
A lower bound ti for Too can be obtained from Eq. H27() 
by assuming that the initial acceleration is maintained 
until Q{ti)/Qoo = 1 is reached. 



C 

Too > Tl = — T. 

a 



(29) 



A better estimate T2 for Too is obtained from Eq. H28|l by 
demanding Q{t2)/Qoo = 1 - e"^. 



Too ~ 'r2 = T. 

ai 



(30) 
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FIG. 1: A log-log plot of the flow rate Q{t)/Qoo as a function 
of time t/r. The dashed line is the short-time approxima- 
tion Eq. (1271 . while the dashed-dotted line is the long-time 
approximation Eq. 12811 . both for the case of a circle, i.e., us- 
ing C/a = 2 and ai = 1.45 as listed in Table |T] The data 
points are the results of time- dependent finite-element simu- 
lations for the cases of the cross section being a circle (white 
circles), a square (gray squares), and an equilateral triangle 
(black triangles) . 



ever, in general the equations have to be solved numer- 
ically, and for this purpose we have used the commer- 
cially available finite-element software Comsol 3.2 (see 
www.comsol.com). Numbers for a selection of geometries 
are tabulated in Table 

The circle is the most compact shape and consequently 
it has the largest value for Ai/A, i.e., the mode has 
the relatively largest spatial occupation of the total area. 
The eigenvalue ai is of the order unity for compact shapes 
and in general it tends to increase slightly with increas- 
ing values of C. The modest variation from geometry to 
geometry in both ai and the other parameters suggests 
that the dynamics of Q{t) will appear almost universal. 

In order to illustrate the validity of our two asymptotic 
expressions, Eqs. H27|l and (|28|l . we have compared them 
using the values for a circular shape to time-dependent 
finite-element simulations of Eq. As illustrated in 

Fig. n we find a perfect agreement between the asymp- 
totic expressions Eqs. (|27|) and (|28|l and the numerically 
exact data for a circle, a square, and an equilateral tri- 
angle. Comparing the corresponding parameters in Ta- 
ble we would expect all data to almost coincide, which 
is indeed also observed in Fig. ^ The small spread in 
eigenvalues and other parameters thus gives rise to close- 
to-universal dynamics. From the plot it is also clear that 
T is indeed a good estimate for the time it takes to reach 
the steady state. 



Using the parameter values for the circle listed in Table |T1 
we find the values 

ri =0.5t<T2 = 2.1t«Too. (31) 



VII. NUMERICAL RESULTS 

Only few geometries allow analytical solutions of both 
the Helmholz equation and the Poisson equation. The 
circle is of course the most well-known example, but 
the equilateral triangle is another exception. How- 



VIII. CONCLUSIONS 

In conclusion, by using a compact Hilbert space for- 
malism we have shown how the initial dynamics in the 
onset of Poiseuille flow is governed by a universal linear 
raise in flow rate Q{t) over a universal time-scale t above 
which it saturates exponentially to the steady-state value 
Qoo- The steady state is reached after a time Too ~ Ct / a. 
Apart from being a fascinating example of universal dy- 
namics for a complex problem our results may have im- 
portant applications in design of real-time programmable 
pressure-driven micro-fluidic networks. 
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